Проект IV

Проект IV#

Создание интерактивной карты#

1. Тема карты#

🛫 Отображение выполненных работ по аэрофотосъемке с помощью Folium 🛬#

Цель:

Создание интерактивной карты, отображающую географию и некоторую статистику проведенных работ по аэрофотосъемке и воздушному лазерному сканированию за временной срез 2016-2024 годы.

Исходя из требований, создам:

4 разных тематических слоя и 2 картографических подложки

Cartodbpositron
OSM
слой Субъекты РФ
слой Работали здесь
слой Места съемок
слой Количество проектов

разных способа изображения данных:

картограмма (choropleth)
катодиаграмма (маркеры с кругами)
кластеризация маркеров
маркеры

другие обязательные элементы:

Layer Control (управление слоями)

Дополнительные интерактивные инструменты:**

Mouse Position
Mini Map
Marker cluster

2. Изучение данных#

2.1 Импорт библиотек#

import folium
import pandas as pd
import geopandas as gpd
import osmnx as ox
import gzip as gz
import json as js


import branca.colormap as cm

from branca.element import Html, MacroElement
from folium.plugins import MousePosition
from folium.plugins import MarkerCluster
from folium.plugins import MiniMap
from folium import IFrame

2.2 Загрузка данных о проектах#

# Загрузим данные о выполненных проектах
df = gpd.read_file('2016-2024_laser_cover.shp')

# Определяем СК
df_crs = df.crs

# Проверяем СК
print("Cистема координат слоя", df_crs)

---------------------------------------------------------------------------
DataSourceError                           Traceback (most recent call last)
Cell In[2], line 2
# Загрузим данные о выполненных проектах
----> 2 df = gpd.read_file('2016-2024_laser_cover.shp')
# Определяем СК
df_crs = df.crs

File /opt/anaconda3/lib/python3.12/site-packages/geopandas/io/file.py:294, in _read_file(filename, bbox, mask, columns, rows, engine, **kwargs)
           from_bytes = True
if engine == "pyogrio":
--> 294     return _read_file_pyogrio(
       filename, bbox=bbox, mask=mask, columns=columns, rows=rows, **kwargs
   )
elif engine == "fiona":
   if pd.api.types.is_file_like(filename):

File /opt/anaconda3/lib/python3.12/site-packages/geopandas/io/file.py:547, in _read_file_pyogrio(path_or_bytes, bbox, mask, rows, **kwargs)
   warnings.warn(
       "The 'include_fields' and 'ignore_fields' keywords are deprecated, and "
       "will be removed in a future release. You can use the 'columns' keyword "
   (...)
       stacklevel=3,
   )
   kwargs["columns"] = kwargs.pop("include_fields")
--> 547 return pyogrio.read_dataframe(path_or_bytes, bbox=bbox, **kwargs)

File /opt/anaconda3/lib/python3.12/site-packages/pyogrio/geopandas.py:265, in read_dataframe(path_or_buffer, layer, encoding, columns, read_geometry, force_2d, skip_features, max_features, where, bbox, mask, fids, sql, sql_dialect, fid_as_index, use_arrow, on_invalid, arrow_to_pandas_kwargs, **kwargs)
if not use_arrow:
   # For arrow, datetimes are read as is.
   # For numpy IO, datetimes are read as string values to preserve timezone info
   # as numpy does not directly support timezones.
   kwargs["datetime_as_string"] = True
--> 265 result = read_func(
   path_or_buffer,
   layer=layer,
   encoding=encoding,
   columns=columns,
   read_geometry=read_geometry,
   force_2d=gdal_force_2d,
   skip_features=skip_features,
   max_features=max_features,
   where=where,
   bbox=bbox,
   mask=mask,
   fids=fids,
   sql=sql,
   sql_dialect=sql_dialect,
   return_fids=fid_as_index,
   **kwargs,
)
if use_arrow:
   meta, table = result

File /opt/anaconda3/lib/python3.12/site-packages/pyogrio/raw.py:198, in read(path_or_buffer, layer, encoding, columns, read_geometry, force_2d, skip_features, max_features, where, bbox, mask, fids, sql, sql_dialect, return_fids, datetime_as_string, **kwargs)
"""Read OGR data source into numpy arrays.

IMPORTANT: non-linear geometry types (e.g., MultiSurface) are converted
   (...)

"""
dataset_kwargs = _preprocess_options_key_value(kwargs) if kwargs else {}
--> 198 return ogr_read(
   get_vsi_path_or_buffer(path_or_buffer),
   layer=layer,
   encoding=encoding,
   columns=columns,
   read_geometry=read_geometry,
   force_2d=force_2d,
   skip_features=skip_features,
   max_features=max_features or 0,
   where=where,
   bbox=bbox,
   mask=_mask_to_wkb(mask),
   fids=fids,
   sql=sql,
   sql_dialect=sql_dialect,
   return_fids=return_fids,
   dataset_kwargs=dataset_kwargs,
   datetime_as_string=datetime_as_string,
)

File pyogrio/_io.pyx:1240, in pyogrio._io.ogr_read()

File pyogrio/_io.pyx:220, in pyogrio._io.ogr_open()

DataSourceError: 2016-2024_laser_cover.shp: No such file or directory

# Посмотрим на таблицу
df.head(3)

Расшифрую наполнение таблицы:

Столбцы	значения
OBJECTID	Числовой идентификатор объекта
project_id	Внутренне название проекта
area	Площадь в км²
year	Год съемки
SK	Итоговая система координат
model	Наличие 3д модели (есть/нет)
object_id	Внутренне название объекта съемки
Shape_Leng	Не обращаем внимания
Shape_Area	Не обращаем внимания
geometry	Координаты в градусах

Для дальнейшей работы мне понадобятся центроиды, по-этому нужно спроецировать градусы в метры. Сложность в том, что данные покрывают всютерриторию России. Остановимся на Web Mercator, для наших целей точности достаточно.

# Перепроецируем в спроецированную систему координат, выберем Web Mercator
df_webmer = df.to_crs(epsg=3857)

# Создаем GeoDataFrame с центроидами
centroids = df_webmer.copy()
centroids['geometry'] = centroids.geometry.centroid  # заменяем геометрию на центроиды

# Отображаем полигоны и центроиды на одной карте
m = df_webmer.explore(color='red', name='Полигоны съемки')  # полигоны
centroids.explore(m=m, 
                  color='blue', 
                  marker_kwds={'radius': 2}, 
                  name='Центроиды'
                  )  # точки

# Покажем слои
folium.LayerControl().add_to(m)
m

centroids.head(3)

# Переведем в систему координат обратно WGS84
centroids_wgs84 = centroids.to_crs(epsg=4326)

2.3 Загрузка данных о субъектах РФ#

Так как изначально файл превышает 25 МБ, он был заархивирован и, затем, преобразован в GeoDataFrame

# Путь к сжатому файлу
input_path = "russia_region.geojson.gz"

# Чтение и преобразование в GeoDataFrame
with gz.open(input_path, "rt", encoding="utf-8") as f: # "rt" - режим чтения текста (r = read, t = text)
    geojson_data = js.load(f)

# Преобразуем в GeoDataFrame
admin = gpd.GeoDataFrame.from_features(geojson_data["features"])

# Указываем систему координат, если она известна (например, WGS84)
admin.set_crs("EPSG:4326", inplace=True)  # ← Замените, если ваша СК другая

# Проверки
print("Система координат слоя:", admin.crs)
print("Тип данных слоя:", admin.geom_type.unique())

# Отобразим данные на карте

# Создадим базовую карту
m = folium.Map(tiles='CartoDB positron')

# Добавим слой на карту

admin.explore(m=m, color='lightblue', opacity=0.8, name='субъекты РФ')  # полигоны
df.explore(m=m, color='red', name='данные съемки')
centroids_wgs84.explore(m=m, color='blue', marker_kwds={'radius': 2}, name='Центроиды')  # точки

# Добавляем контроль слоев для возможности включения/выключения
folium.LayerControl().add_to(m)

# Автоматически подстраиваем масштаб под все объекты
m.fit_bounds(m.get_bounds())

# Отображаем карту
m

3. Визуализация данных#

Мы создадим хороплетную карту, где будет показано сколько проектов снято в регионе РФ. Для этого нам понадобится сделать пространственный join центроидов и регионов.

# Делаем пространственный join: для каждой точки подбираем тот регион , в который она попадает
#    predicate='within' означает, что точка должна лежать внутри полигона округа
join_map = gpd.sjoin(centroids_wgs84, admin, how='left', predicate='within')

# Переименовываем столбцы, чтобы с ними было удобнее работать 
join_map = join_map.rename(columns={
    'name': 'region_name',        # новое имя для столбца "название округа"
})

join_map.head(5)

join_map.columns

project_counts = join_map.groupby('region_name').agg(
    project_count=('region_name', 'size'),
    total_area=('area', 'sum')  # округляем до двух знаков после запятой
).reset_index()
project_counts['total_area'] = project_counts['total_area'].round(0).astype(int)  # округляем до целого числа
project_counts

# # Группировка по region_name и подсчет количества проектов
# project_counts = join_map.groupby('region_name').size().reset_index(name='project_count')
# project_counts

Так как при пространственном объединении (sjoin), геометрия берется от левой таблицы (centroids_wgs84), то есть — точки и join_map содержит геометрию точек, а не полигонов, нужно снова присоединить информацию о полигонах (регионов) к аггрегированной таблице project_counts чтобы построить хороплетную карту.

# Проверяем тип данных
print("Тип данных слоя", join_map.geom_type.unique())

# Присоединяем геометрию регионов обратно (из admin)
choropleth_data = admin.merge(project_counts, left_on='name', right_on='region_name', how='left')
# Проверяем тип данных
print("Тип данных слоя", choropleth_data.geom_type.unique())

# Если у каких-то регионов нет точек в project_count будет NaN.
# Заменим нулями и переведем в целочисленный тип
choropleth_data['project_count'] = choropleth_data['project_count'].fillna(0).astype(int)
choropleth_data['project_count'].isna().sum()

# Удалим лишние столбцы
choropleth_data = choropleth_data.drop(columns=['NUMPOINTS', 'region_name'])
# Заменим NaN значения в столбце total_area на 0
choropleth_data['total_area'] = choropleth_data['total_area'].fillna(0).astype(int)

choropleth_data.head()

Отметим, что для визуализации будем использовать name, потому что там полные данные

3.1 Создаем хороплетную карту#

# Создаем карту с пустым фоном
m = folium.Map(location=[68.13, 98.43], 
               zoom_start=2.5, 
               tiles=None )

# Добавляем подложки как отключаемые слои
folium.TileLayer(
    tiles='cartodbpositron',
    name='Cartodbpositron',
    control=True
).add_to(m)

folium.TileLayer(
    tiles='OpenStreetMap',
    name='OpenStreetMap',
    control=True
).add_to(m)

# Добавляем слой с субъектами
folium.GeoJson(
    choropleth_data,
    name='Субьекты РФ',
    style_function=lambda feature: {
        'fillColor': '#4CAFE5',
        'color': '#FAFAFA',
        'weight': 1.5,
        'fillOpacity': 0.4
    },
    tooltip=folium.GeoJsonTooltip(
        fields=['name'],
        aliases=[''],
        localize=True,
        sticky=True,
        style="padding-left: 4px; text-align: left;"
    )
).add_to(m)

# Добавляем слой: работали здесь (только > 0)
folium.GeoJson(
    choropleth_data[choropleth_data['project_count'] > 0],
    name='Работали здесь',
    style_function=lambda feature: {
        'fillColor': '#013A63',
        'color': '#FAFAFA',
        'weight': 1.5,
        'fillOpacity': 0.7
    }
).add_to(m)

# Создаем FeatureGroup для слоя с маркерами (изначально выключен)
marker_layer = folium.FeatureGroup(name='Места съемок', show=False)

# Кластер для маркеров
marker_cluster = MarkerCluster().add_to(marker_layer)

# Добавляем точки в кластер
for idx, row in centroids_wgs84.iterrows():
    location = [row.geometry.y, row.geometry.x]
    popup_text = f"Год съемки: {row['year']}"
    
    folium.Marker(
        location=location,
        popup=popup_text,
        icon=folium.Icon(color='blue', icon='plane', prefix='fa')
    ).add_to(marker_cluster)

# Добавляем слой с маркерами на карту
marker_layer.add_to(m)

<folium.map.FeatureGroup at 0x169f3f0b0>

# Считаем центроиды (НЕ сохраняем в geometry!)
centroids = choropleth_data.geometry.centroid
choropleth_data['lon'] = centroids.x
choropleth_data['lat'] = centroids.y

# Создаем пустой FeatureGroup для кругов (по умолчанию скрыт)
circle_layer = folium.FeatureGroup(name='Количество проектов')

# Добавляем круги
for _, row in choropleth_data.iterrows():
    if row['project_count'] > 0:
        html = f"""
        <div style="text-align: left; font-weight: normal; font-size: 12px;">
            <strong>{row['name']}</strong><br>
            Проектов: {row['project_count']}<br>
            Площадь съемки: {row['total_area']} км²
        </div>
        """
        iframe = IFrame(html=html, width=200, height=70)
        popup = folium.Popup(iframe, max_width=200)

        folium.CircleMarker(
            location=[row['lat'], row['lon']],
            radius=row['project_count'] ** 0.5 * 2,
            color='#01497C',
            weight=1,
            fill=True,
            fill_color="#CD8F0B",
            fill_opacity=0.8,
            popup=popup
        ).add_to(circle_layer)

circle_layer.add_to(m)

# Добавляем мини-карту
MiniMap().add_to(m)

# Добавляем координаты и котроль слоев
MousePosition().add_to(m)

folium.LayerControl(collapsed=False).add_to(m)

# Добавим заголовок к карте
title_html = '''
     <h3 align="center" style="font-size:17px">
     <b>Выполненные работы по аэрофотосъемке и воздушному лазерному сканированию в субъектах РФ</b></h3>
     '''
m.get_root().html.add_child(folium.Element(title_html))

<branca.element.Element at 0x1253e2b10>

Make this Notebook Trusted to load map: File -> Trust Notebook

m.save("index.html")